# ==================== 1. 加载包 ====================
library(caret)      # 用于 createFolds 和 confusionMatrix
library(pROC)
library(ggplot2)

# ==================== 2. 使用您已读入的数据 ====================
train_data <- Ttrain
test_data  <- Ttest

cat("训练集样本量：", nrow(train_data), "\n")
cat("测试集样本量：", nrow(test_data), "\n")

# ==================== 3. 数据预处理 ====================

# 3.1 将结局变量转为因子，水平命名为 "No" 和 "Yes"（阳性为 "Yes"）
train_data$营养不良风险 <- factor(train_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))
test_data$营养不良风险  <- factor(test_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))

# 3.2 将所有字符型变量转为因子
for (col in names(train_data)) {
  if (is.character(train_data[[col]])) {
    train_data[[col]] <- as.factor(train_data[[col]])
    test_data[[col]]  <- as.factor(test_data[[col]])
  }
}

# 3.3 缺失值填补（用训练集的中位数/众数填补）
impute_vals <- list()
for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    impute_vals[[col]] <- median(train_data[[col]], na.rm = TRUE)
  } else if (is.factor(train_data[[col]])) {
    freq <- table(train_data[[col]])
    impute_vals[[col]] <- names(freq)[which.max(freq)]
  }
}

for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  } else if (is.factor(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  }
}

cat("填补后训练集缺失值总数：", sum(is.na(train_data)), "\n")
cat("填补后测试集缺失值总数：", sum(is.na(test_data)), "\n")

# 检查类别分布
cat("\n训练集类别分布：\n")
print(table(train_data$营养不良风险))
cat("测试集类别分布：\n")
print(table(test_data$营养不良风险))

if (length(unique(train_data$营养不良风险)) < 2) {
  stop("错误：训练集中只有一个类别，请检查数据！")
}

# ==================== 4. 手动 10 折交叉验证 ====================
set.seed(123)
folds <- createFolds(train_data$营养不良风险, k = 10, list = TRUE)

# 初始化存储数据框
cv_metrics <- data.frame(
  Fold = 1:10,
  AUC = NA,
  Accuracy = NA,
  Sensitivity = NA,
  Specificity = NA
)

for (i in 1:10) {
  val_idx <- folds[[i]]
  train_fold <- train_data[-val_idx, ]
  val_fold   <- train_data[val_idx, ]
  
  # 训练模型
  model_fold <- glm(营养不良风险 ~ ., data = train_fold, family = binomial)
  
  # 预测概率
  pred_prob <- predict(model_fold, newdata = val_fold, type = "response")
  pred_class <- ifelse(pred_prob > 0.5, "Yes", "No")
  pred_class <- factor(pred_class, levels = c("No", "Yes"))
  
  # 计算混淆矩阵指标
  cm <- confusionMatrix(pred_class, val_fold$营养不良风险, positive = "Yes")
  cv_metrics$Accuracy[i]   <- cm$overall["Accuracy"]
  cv_metrics$Sensitivity[i] <- cm$byClass["Sensitivity"]
  cv_metrics$Specificity[i] <- cm$byClass["Specificity"]
  
  # 计算 AUC
  roc_obj <- roc(response = val_fold$营养不良风险, predictor = pred_prob, 
                 levels = c("No", "Yes"))
  cv_metrics$AUC[i] <- auc(roc_obj)
}

# 输出 CV 结果
cat("\n========== 10 折交叉验证每折结果 ==========\n")
print(cv_metrics)

cat("\n========== 10 折交叉验证汇总（均值 ± 标准差） ==========\n")
cat(sprintf("AUC          : %.4f ± %.4f\n", mean(cv_metrics$AUC), sd(cv_metrics$AUC)))
cat(sprintf("准确率 (Acc) : %.4f ± %.4f\n", mean(cv_metrics$Accuracy), sd(cv_metrics$Accuracy)))
cat(sprintf("敏感度 (Sens): %.4f ± %.4f\n", mean(cv_metrics$Sensitivity), sd(cv_metrics$Sensitivity)))
cat(sprintf("特异度 (Spec): %.4f ± %.4f\n", mean(cv_metrics$Specificity), sd(cv_metrics$Specificity)))

# ==================== 5. 用全量训练集训练最终模型 ====================
set.seed(123)
lr_model <- glm(营养不良风险 ~ ., data = train_data, family = binomial)

# ==================== 6. 测试集预测 ====================
pred_prob_pos <- predict(lr_model, newdata = test_data, type = "response")
pred_class <- ifelse(pred_prob_pos > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = c("No", "Yes"))

# ==================== 7. 保存测试集预测结果（供合并 ROC） ====================
test_results <- data.frame(
  true_label = test_data$营养不良风险,
  pred_prob  = pred_prob_pos
)
write.csv(test_results, file = "lr_test_predictions.csv", row.names = FALSE)
cat("\n✅ 测试集预测结果已保存为：lr_test_predictions.csv\n")

# ==================== 8. 测试集完整评价指标 ====================
cat("\n========== 测试集模型性能（完整指标） ==========\n")

cm <- confusionMatrix(pred_class, test_data$营养不良风险, positive = "Yes")
cat("混淆矩阵：\n")
print(cm$table)

cat("\n--- 总体指标 ---\n")
cat(sprintf("准确率 (Accuracy)          : %.2f%%  (95%% CI: %.2f%% - %.2f%%)\n", 
            cm$overall["Accuracy"] * 100,
            cm$overall["AccuracyLower"] * 100,
            cm$overall["AccuracyUpper"] * 100))
cat(sprintf("Kappa 值                  : %.4f\n", cm$overall["Kappa"]))

cat("\n--- 类别层面指标（以阳性类别 \"Yes\" 为关注） ---\n")
cat(sprintf("敏感度 (Sensitivity/Recall): %.2f%%\n", cm$byClass["Sensitivity"] * 100))
cat(sprintf("特异度 (Specificity)       : %.2f%%\n", cm$byClass["Specificity"] * 100))
cat(sprintf("阳性预测值 (Precision)     : %.2f%%\n", cm$byClass["Pos Pred Value"] * 100))
cat(sprintf("阴性预测值 (Neg Pred Value): %.2f%%\n", cm$byClass["Neg Pred Value"] * 100))
cat(sprintf("F1 分数                   : %.4f\n", cm$byClass["F1"]))
cat(sprintf("平衡准确率 (Balanced Acc)  : %.2f%%\n", cm$byClass["Balanced Accuracy"] * 100))

# ROC 和 AUC
roc_obj <- roc(response = test_data$营养不良风险, predictor = pred_prob_pos, 
               levels = c("No", "Yes"))
auc_val <- auc(roc_obj)
ci_val <- ci(roc_obj)
cat(sprintf("\nAUC：%.3f (95%% CI: %.3f - %.3f)\n", auc_val, ci_val[1], ci_val[3]))

# ==================== 9. 绘制并保存 ROC 曲线 ====================
roc_df <- data.frame(
  FPR = 1 - roc_obj$specificities,
  TPR = roc_obj$sensitivities
)

p <- ggplot(roc_df, aes(x = FPR, y = TPR)) +
  geom_line(color = "#2E86C1", size = 1.2) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  labs(title = "逻辑回归 ROC 曲线 (测试集)",
       subtitle = paste0("AUC = ", round(auc_val, 3), 
                         " (95% CI: ", round(ci_val[1], 3), " - ", round(ci_val[3], 3), ")"),
       x = "假阳性率 (1 - 特异度)",
       y = "真阳性率 (敏感度)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10)) +
  coord_equal() +
  xlim(0, 1) + ylim(0, 1)

ggsave("lr_roc_curve.pdf", plot = p, width = 6, height = 5, dpi = 300)
ggsave("lr_roc_curve.png", plot = p, width = 6, height = 5, dpi = 300)
print(p)

cat("\n✅ ROC 曲线已保存为：lr_roc_curve.pdf 和 lr_roc_curve.png\n")
cat("\n✅ 逻辑回归建模完成（手动10折CV + 完整指标 + ROC图）！\n")




library(caret)      # createFolds, confusionMatrix
library(e1071)      # svm
library(pROC)
library(ggplot2)

# ==================== 2. 使用您已读入的数据 ====================
train_data <- Ttrain
test_data  <- Ttest

cat("训练集样本量：", nrow(train_data), "\n")
cat("测试集样本量：", nrow(test_data), "\n")

# ==================== 3. 数据预处理 ====================

# 3.1 将结局变量转为因子，水平命名为 "No" 和 "Yes"（阳性为 "Yes"）
train_data$营养不良风险 <- factor(train_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))
test_data$营养不良风险  <- factor(test_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))

# 3.2 将所有字符型变量转为因子（SVM 中可直接使用因子）
for (col in names(train_data)) {
  if (is.character(train_data[[col]])) {
    train_data[[col]] <- as.factor(train_data[[col]])
    test_data[[col]]  <- as.factor(test_data[[col]])
  }
}

# 3.3 缺失值填补（与逻辑回归一致）
impute_vals <- list()
for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    impute_vals[[col]] <- median(train_data[[col]], na.rm = TRUE)
  } else if (is.factor(train_data[[col]])) {
    freq <- table(train_data[[col]])
    impute_vals[[col]] <- names(freq)[which.max(freq)]
  }
}

for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  } else if (is.factor(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  }
}

cat("填补后训练集缺失值总数：", sum(is.na(train_data)), "\n")
cat("填补后测试集缺失值总数：", sum(is.na(test_data)), "\n")

cat("\n训练集类别分布：\n")
print(table(train_data$营养不良风险))
cat("测试集类别分布：\n")
print(table(test_data$营养不良风险))

if (length(unique(train_data$营养不良风险)) < 2) {
  stop("错误：训练集中只有一个类别，无法建模")
}

# ==================== 4. 手动 10 折交叉验证 ====================
set.seed(123)
folds <- createFolds(train_data$营养不良风险, k = 10, list = TRUE)

cv_metrics <- data.frame(
  Fold = 1:10,
  AUC = NA,
  Accuracy = NA,
  Sensitivity = NA,
  Specificity = NA
)

for (i in 1:10) {
  val_idx <- folds[[i]]
  train_fold <- train_data[-val_idx, ]
  val_fold   <- train_data[val_idx, ]
  
  # 训练 SVM（使用线性核，可改为 "radial" 等）
  model_fold <- svm(营养不良风险 ~ ., 
                    data = train_fold, 
                    kernel = "linear", 
                    probability = TRUE,
                    scale = TRUE)
  
  # 预测概率
  pred_prob_obj <- predict(model_fold, newdata = val_fold, probability = TRUE)
  pred_prob <- attr(pred_prob_obj, "probabilities")[, "Yes"]   # 提取 "Yes" 的概率
  
  pred_class <- ifelse(pred_prob > 0.5, "Yes", "No")
  pred_class <- factor(pred_class, levels = c("No", "Yes"))
  
  # 计算指标
  cm <- confusionMatrix(pred_class, val_fold$营养不良风险, positive = "Yes")
  cv_metrics$Accuracy[i]   <- cm$overall["Accuracy"]
  cv_metrics$Sensitivity[i] <- cm$byClass["Sensitivity"]
  cv_metrics$Specificity[i] <- cm$byClass["Specificity"]
  
  roc_obj <- roc(response = val_fold$营养不良风险, predictor = pred_prob, 
                 levels = c("No", "Yes"))
  cv_metrics$AUC[i] <- auc(roc_obj)
}

cat("\n========== 10 折交叉验证每折结果 ==========\n")
print(cv_metrics)

cat("\n========== 10 折交叉验证汇总（均值 ± 标准差） ==========\n")
cat(sprintf("AUC          : %.4f ± %.4f\n", mean(cv_metrics$AUC), sd(cv_metrics$AUC)))
cat(sprintf("准确率 (Acc) : %.4f ± %.4f\n", mean(cv_metrics$Accuracy), sd(cv_metrics$Accuracy)))
cat(sprintf("敏感度 (Sens): %.4f ± %.4f\n", mean(cv_metrics$Sensitivity), sd(cv_metrics$Sensitivity)))
cat(sprintf("特异度 (Spec): %.4f ± %.4f\n", mean(cv_metrics$Specificity), sd(cv_metrics$Specificity)))

# ==================== 5. 用全量训练集训练最终 SVM 模型 ====================
set.seed(123)
svm_model <- svm(营养不良风险 ~ ., 
                 data = train_data, 
                 kernel = "linear", 
                 probability = TRUE,
                 scale = TRUE)

# ==================== 6. 测试集预测 ====================
pred_prob_obj <- predict(svm_model, newdata = test_data, probability = TRUE)
pred_prob_pos <- attr(pred_prob_obj, "probabilities")[, "Yes"]

pred_class <- ifelse(pred_prob_pos > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = c("No", "Yes"))

# ==================== 7. 保存测试集预测结果 ====================
test_results <- data.frame(
  true_label = test_data$营养不良风险,
  pred_prob  = pred_prob_pos
)
write.csv(test_results, file = "svm_test_predictions.csv", row.names = FALSE)
cat("\n✅ 测试集预测结果已保存为：svm_test_predictions.csv\n")

# ==================== 8. 测试集完整评价指标 ====================
cat("\n========== 测试集模型性能（完整指标） ==========\n")

cm <- confusionMatrix(pred_class, test_data$营养不良风险, positive = "Yes")
cat("混淆矩阵：\n")
print(cm$table)

cat("\n--- 总体指标 ---\n")
cat(sprintf("准确率 (Accuracy)          : %.2f%%  (95%% CI: %.2f%% - %.2f%%)\n", 
            cm$overall["Accuracy"] * 100,
            cm$overall["AccuracyLower"] * 100,
            cm$overall["AccuracyUpper"] * 100))
cat(sprintf("Kappa 值                  : %.4f\n", cm$overall["Kappa"]))

cat("\n--- 类别层面指标（以阳性类别 \"Yes\" 为关注） ---\n")
cat(sprintf("敏感度 (Sensitivity/Recall): %.2f%%\n", cm$byClass["Sensitivity"] * 100))
cat(sprintf("特异度 (Specificity)       : %.2f%%\n", cm$byClass["Specificity"] * 100))
cat(sprintf("阳性预测值 (Precision)     : %.2f%%\n", cm$byClass["Pos Pred Value"] * 100))
cat(sprintf("阴性预测值 (Neg Pred Value): %.2f%%\n", cm$byClass["Neg Pred Value"] * 100))
cat(sprintf("F1 分数                   : %.4f\n", cm$byClass["F1"]))
cat(sprintf("平衡准确率 (Balanced Acc)  : %.2f%%\n", cm$byClass["Balanced Accuracy"] * 100))

roc_obj <- roc(response = test_data$营养不良风险, predictor = pred_prob_pos, 
               levels = c("No", "Yes"))
auc_val <- auc(roc_obj)
ci_val <- ci(roc_obj)
cat(sprintf("\nAUC：%.3f (95%% CI: %.3f - %.3f)\n", auc_val, ci_val[1], ci_val[3]))

# ==================== 9. 绘制并保存 ROC 曲线 ====================
roc_df <- data.frame(
  FPR = 1 - roc_obj$specificities,
  TPR = roc_obj$sensitivities
)

p <- ggplot(roc_df, aes(x = FPR, y = TPR)) +
  geom_line(color = "#E74C3C", size = 1.2) +   # 用红色区分逻辑回归的蓝色
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  labs(title = "SVM ROC 曲线 (测试集)",
       subtitle = paste0("AUC = ", round(auc_val, 3), 
                         " (95% CI: ", round(ci_val[1], 3), " - ", round(ci_val[3], 3), ")"),
       x = "假阳性率 (1 - 特异度)",
       y = "真阳性率 (敏感度)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10)) +
  coord_equal() +
  xlim(0, 1) + ylim(0, 1)

ggsave("svm_roc_curve.pdf", plot = p, width = 6, height = 5, dpi = 300)
ggsave("svm_roc_curve.png", plot = p, width = 6, height = 5, dpi = 300)
print(p)

cat("\n✅ ROC 曲线已保存为：svm_roc_curve.pdf 和 svm_roc_curve.png\n")
cat("\n✅ SVM 建模完成（手动10折CV + 完整指标 + ROC图）！\n")


# ==================== 1. 加载包 ====================
library(caret)      # createFolds, confusionMatrix
library(rpart)      # 决策树
library(pROC)
library(ggplot2)

# ==================== 2. 使用您已读入的数据 ====================
train_data <- Ttrain
test_data  <- Ttest

cat("训练集样本量：", nrow(train_data), "\n")
cat("测试集样本量：", nrow(test_data), "\n")

# ==================== 3. 数据预处理 ====================

# 3.1 将结局变量转为因子，水平命名为 "No" 和 "Yes"（阳性为 "Yes"）
train_data$营养不良风险 <- factor(train_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))
test_data$营养不良风险  <- factor(test_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))

# 3.2 将所有字符型变量转为因子（rpart 能处理因子，但需确保）
for (col in names(train_data)) {
  if (is.character(train_data[[col]])) {
    train_data[[col]] <- as.factor(train_data[[col]])
    test_data[[col]]  <- as.factor(test_data[[col]])
  }
}

# 3.3 缺失值填补（与之前一致，用训练集的中位数/众数填补）
impute_vals <- list()
for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    impute_vals[[col]] <- median(train_data[[col]], na.rm = TRUE)
  } else if (is.factor(train_data[[col]])) {
    freq <- table(train_data[[col]])
    impute_vals[[col]] <- names(freq)[which.max(freq)]
  }
}

for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  } else if (is.factor(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  }
}

cat("填补后训练集缺失值总数：", sum(is.na(train_data)), "\n")
cat("填补后测试集缺失值总数：", sum(is.na(test_data)), "\n")

cat("\n训练集类别分布：\n")
print(table(train_data$营养不良风险))
cat("测试集类别分布：\n")
print(table(test_data$营养不良风险))

if (length(unique(train_data$营养不良风险)) < 2) {
  stop("错误：训练集中只有一个类别，无法建模")
}

# ==================== 4. 手动 10 折交叉验证 ====================
set.seed(123)
folds <- createFolds(train_data$营养不良风险, k = 10, list = TRUE)

cv_metrics <- data.frame(
  Fold = 1:10,
  AUC = NA,
  Accuracy = NA,
  Sensitivity = NA,
  Specificity = NA
)

for (i in 1:10) {
  val_idx <- folds[[i]]
  train_fold <- train_data[-val_idx, ]
  val_fold   <- train_data[val_idx, ]
  
  # 训练决策树（使用默认参数，可添加 control = rpart.control(minsplit = 20, cp = 0.01)）
  model_fold <- rpart(营养不良风险 ~ ., 
                      data = train_fold, 
                      method = "class")   # 分类树
  
  # 预测概率（返回矩阵，列名为类别）
  pred_prob_mat <- predict(model_fold, newdata = val_fold, type = "prob")
  pred_prob <- pred_prob_mat[, "Yes"]   # 提取阳性类别概率
  
  pred_class <- ifelse(pred_prob > 0.5, "Yes", "No")
  pred_class <- factor(pred_class, levels = c("No", "Yes"))
  
  # 计算指标
  cm <- confusionMatrix(pred_class, val_fold$营养不良风险, positive = "Yes")
  cv_metrics$Accuracy[i]   <- cm$overall["Accuracy"]
  cv_metrics$Sensitivity[i] <- cm$byClass["Sensitivity"]
  cv_metrics$Specificity[i] <- cm$byClass["Specificity"]
  
  roc_obj <- roc(response = val_fold$营养不良风险, predictor = pred_prob, 
                 levels = c("No", "Yes"))
  cv_metrics$AUC[i] <- auc(roc_obj)
}

cat("\n========== 10 折交叉验证每折结果 ==========\n")
print(cv_metrics)

cat("\n========== 10 折交叉验证汇总（均值 ± 标准差） ==========\n")
cat(sprintf("AUC          : %.4f ± %.4f\n", mean(cv_metrics$AUC), sd(cv_metrics$AUC)))
cat(sprintf("准确率 (Acc) : %.4f ± %.4f\n", mean(cv_metrics$Accuracy), sd(cv_metrics$Accuracy)))
cat(sprintf("敏感度 (Sens): %.4f ± %.4f\n", mean(cv_metrics$Sensitivity), sd(cv_metrics$Sensitivity)))
cat(sprintf("特异度 (Spec): %.4f ± %.4f\n", mean(cv_metrics$Specificity), sd(cv_metrics$Specificity)))

# ==================== 5. 用全量训练集训练最终决策树模型 ====================
set.seed(123)
dt_model <- rpart(营养不良风险 ~ ., data = train_data, method = "class")

# ==================== 6. 测试集预测 ====================
pred_prob_mat <- predict(dt_model, newdata = test_data, type = "prob")
pred_prob_pos <- pred_prob_mat[, "Yes"]

pred_class <- ifelse(pred_prob_pos > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = c("No", "Yes"))

# ==================== 7. 保存测试集预测结果 ====================
test_results <- data.frame(
  true_label = test_data$营养不良风险,
  pred_prob  = pred_prob_pos
)
write.csv(test_results, file = "dt_test_predictions.csv", row.names = FALSE)
cat("\n✅ 测试集预测结果已保存为：dt_test_predictions.csv\n")

# ==================== 8. 测试集完整评价指标 ====================
cat("\n========== 测试集模型性能（完整指标） ==========\n")

cm <- confusionMatrix(pred_class, test_data$营养不良风险, positive = "Yes")
cat("混淆矩阵：\n")
print(cm$table)

cat("\n--- 总体指标 ---\n")
cat(sprintf("准确率 (Accuracy)          : %.2f%%  (95%% CI: %.2f%% - %.2f%%)\n", 
            cm$overall["Accuracy"] * 100,
            cm$overall["AccuracyLower"] * 100,
            cm$overall["AccuracyUpper"] * 100))
cat(sprintf("Kappa 值                  : %.4f\n", cm$overall["Kappa"]))

cat("\n--- 类别层面指标（以阳性类别 \"Yes\" 为关注） ---\n")
cat(sprintf("敏感度 (Sensitivity/Recall): %.2f%%\n", cm$byClass["Sensitivity"] * 100))
cat(sprintf("特异度 (Specificity)       : %.2f%%\n", cm$byClass["Specificity"] * 100))
cat(sprintf("阳性预测值 (Precision)     : %.2f%%\n", cm$byClass["Pos Pred Value"] * 100))
cat(sprintf("阴性预测值 (Neg Pred Value): %.2f%%\n", cm$byClass["Neg Pred Value"] * 100))
cat(sprintf("F1 分数                   : %.4f\n", cm$byClass["F1"]))
cat(sprintf("平衡准确率 (Balanced Acc)  : %.2f%%\n", cm$byClass["Balanced Accuracy"] * 100))

roc_obj <- roc(response = test_data$营养不良风险, predictor = pred_prob_pos, 
               levels = c("No", "Yes"))
auc_val <- auc(roc_obj)
ci_val <- ci(roc_obj)
cat(sprintf("\nAUC：%.3f (95%% CI: %.3f - %.3f)\n", auc_val, ci_val[1], ci_val[3]))

# ==================== 9. 绘制并保存 ROC 曲线 ====================
roc_df <- data.frame(
  FPR = 1 - roc_obj$specificities,
  TPR = roc_obj$sensitivities
)

p <- ggplot(roc_df, aes(x = FPR, y = TPR)) +
  geom_line(color = "#2CA02C", size = 1.2) +   # 用绿色区分
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  labs(title = "决策树 ROC 曲线 (测试集)",
       subtitle = paste0("AUC = ", round(auc_val, 3), 
                         " (95% CI: ", round(ci_val[1], 3), " - ", round(ci_val[3], 3), ")"),
       x = "假阳性率 (1 - 特异度)",
       y = "真阳性率 (敏感度)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10)) +
  coord_equal() +
  xlim(0, 1) + ylim(0, 1)

ggsave("dt_roc_curve.pdf", plot = p, width = 6, height = 5, dpi = 300)
ggsave("dt_roc_curve.png", plot = p, width = 6, height = 5, dpi = 300)
print(p)

cat("\n✅ ROC 曲线已保存为：dt_roc_curve.pdf 和 dt_roc_curve.png\n")
cat("\n✅ 决策树建模完成（手动10折CV + 完整指标 + ROC图）！\n")


library(caret)          # createFolds, confusionMatrix
library(randomForest)   # 随机森林
library(pROC)
library(ggplot2)

# ==================== 2. 使用您已读入的数据 ====================
train_data <- Ttrain
test_data  <- Ttest

cat("训练集样本量：", nrow(train_data), "\n")
cat("测试集样本量：", nrow(test_data), "\n")

# ==================== 3. 数据预处理 ====================

# 3.1 将结局变量转为因子，水平命名为 "No" 和 "Yes"（阳性为 "Yes"）
train_data$营养不良风险 <- factor(train_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))
test_data$营养不良风险  <- factor(test_data$营养不良风险, 
                                   levels = c(0, 1), 
                                   labels = c("No", "Yes"))

# 3.2 将所有字符型变量转为因子（随机森林需要因子）
for (col in names(train_data)) {
  if (is.character(train_data[[col]])) {
    train_data[[col]] <- as.factor(train_data[[col]])
    test_data[[col]]  <- as.factor(test_data[[col]])
  }
}

# 3.3 缺失值填补（用训练集的中位数/众数填补）
impute_vals <- list()
for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    impute_vals[[col]] <- median(train_data[[col]], na.rm = TRUE)
  } else if (is.factor(train_data[[col]])) {
    freq <- table(train_data[[col]])
    impute_vals[[col]] <- names(freq)[which.max(freq)]
  }
}

for (col in names(train_data)) {
  if (col == "营养不良风险") next
  if (is.numeric(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  } else if (is.factor(train_data[[col]])) {
    val <- impute_vals[[col]]
    train_data[[col]][is.na(train_data[[col]])] <- val
    test_data[[col]][is.na(test_data[[col]])] <- val
  }
}

cat("填补后训练集缺失值总数：", sum(is.na(train_data)), "\n")
cat("填补后测试集缺失值总数：", sum(is.na(test_data)), "\n")

cat("\n训练集类别分布：\n")
print(table(train_data$营养不良风险))
cat("测试集类别分布：\n")
print(table(test_data$营养不良风险))

if (length(unique(train_data$营养不良风险)) < 2) {
  stop("错误：训练集中只有一个类别，无法建模")
}

# ==================== 4. 手动 10 折交叉验证 ====================
set.seed(123)
folds <- createFolds(train_data$营养不良风险, k = 10, list = TRUE)

cv_metrics <- data.frame(
  Fold = 1:10,
  AUC = NA,
  Accuracy = NA,
  Sensitivity = NA,
  Specificity = NA
)

for (i in 1:10) {
  val_idx <- folds[[i]]
  train_fold <- train_data[-val_idx, ]
  val_fold   <- train_data[val_idx, ]
  
  # 训练随机森林（ntree=500 是常用值，可根据运行时间调整）
  model_fold <- randomForest(营养不良风险 ~ ., 
                             data = train_fold, 
                             ntree = 500,
                             importance = FALSE,   # 提升速度
                             mtry = sqrt(ncol(train_fold) - 1))   # 常用默认值
  
  # 预测概率（返回矩阵，列名为类别）
  pred_prob_mat <- predict(model_fold, newdata = val_fold, type = "prob")
  pred_prob <- pred_prob_mat[, "Yes"]   # 提取阳性类别概率
  
  pred_class <- ifelse(pred_prob > 0.5, "Yes", "No")
  pred_class <- factor(pred_class, levels = c("No", "Yes"))
  
  # 计算指标
  cm <- confusionMatrix(pred_class, val_fold$营养不良风险, positive = "Yes")
  cv_metrics$Accuracy[i]   <- cm$overall["Accuracy"]
  cv_metrics$Sensitivity[i] <- cm$byClass["Sensitivity"]
  cv_metrics$Specificity[i] <- cm$byClass["Specificity"]
  
  roc_obj <- roc(response = val_fold$营养不良风险, predictor = pred_prob, 
                 levels = c("No", "Yes"))
  cv_metrics$AUC[i] <- auc(roc_obj)
}

cat("\n========== 10 折交叉验证每折结果 ==========\n")
print(cv_metrics)

cat("\n========== 10 折交叉验证汇总（均值 ± 标准差） ==========\n")
cat(sprintf("AUC          : %.4f ± %.4f\n", mean(cv_metrics$AUC), sd(cv_metrics$AUC)))
cat(sprintf("准确率 (Acc) : %.4f ± %.4f\n", mean(cv_metrics$Accuracy), sd(cv_metrics$Accuracy)))
cat(sprintf("敏感度 (Sens): %.4f ± %.4f\n", mean(cv_metrics$Sensitivity), sd(cv_metrics$Sensitivity)))
cat(sprintf("特异度 (Spec): %.4f ± %.4f\n", mean(cv_metrics$Specificity), sd(cv_metrics$Specificity)))

# ==================== 5. 用全量训练集训练最终随机森林模型 ====================
set.seed(123)
rf_model <- randomForest(营养不良风险 ~ ., 
                         data = train_data, 
                         ntree = 500,
                         importance = FALSE,
                         mtry = sqrt(ncol(train_data) - 1))

# ==================== 6. 测试集预测 ====================
pred_prob_mat <- predict(rf_model, newdata = test_data, type = "prob")
pred_prob_pos <- pred_prob_mat[, "Yes"]

pred_class <- ifelse(pred_prob_pos > 0.5, "Yes", "No")
pred_class <- factor(pred_class, levels = c("No", "Yes"))

# ==================== 7. 保存测试集预测结果 ====================
test_results <- data.frame(
  true_label = test_data$营养不良风险,
  pred_prob  = pred_prob_pos
)
write.csv(test_results, file = "rf_test_predictions.csv", row.names = FALSE)
cat("\n✅ 测试集预测结果已保存为：rf_test_predictions.csv\n")

# ==================== 8. 测试集完整评价指标 ====================
cat("\n========== 测试集模型性能（完整指标） ==========\n")

cm <- confusionMatrix(pred_class, test_data$营养不良风险, positive = "Yes")
cat("混淆矩阵：\n")
print(cm$table)

cat("\n--- 总体指标 ---\n")
cat(sprintf("准确率 (Accuracy)          : %.2f%%  (95%% CI: %.2f%% - %.2f%%)\n", 
            cm$overall["Accuracy"] * 100,
            cm$overall["AccuracyLower"] * 100,
            cm$overall["AccuracyUpper"] * 100))
cat(sprintf("Kappa 值                  : %.4f\n", cm$overall["Kappa"]))

cat("\n--- 类别层面指标（以阳性类别 \"Yes\" 为关注） ---\n")
cat(sprintf("敏感度 (Sensitivity/Recall): %.2f%%\n", cm$byClass["Sensitivity"] * 100))
cat(sprintf("特异度 (Specificity)       : %.2f%%\n", cm$byClass["Specificity"] * 100))
cat(sprintf("阳性预测值 (Precision)     : %.2f%%\n", cm$byClass["Pos Pred Value"] * 100))
cat(sprintf("阴性预测值 (Neg Pred Value): %.2f%%\n", cm$byClass["Neg Pred Value"] * 100))
cat(sprintf("F1 分数                   : %.4f\n", cm$byClass["F1"]))
cat(sprintf("平衡准确率 (Balanced Acc)  : %.2f%%\n", cm$byClass["Balanced Accuracy"] * 100))

roc_obj <- roc(response = test_data$营养不良风险, predictor = pred_prob_pos, 
               levels = c("No", "Yes"))
auc_val <- auc(roc_obj)
ci_val <- ci(roc_obj)
cat(sprintf("\nAUC：%.3f (95%% CI: %.3f - %.3f)\n", auc_val, ci_val[1], ci_val[3]))

# ==================== 9. 绘制并保存 ROC 曲线 ====================
roc_df <- data.frame(
  FPR = 1 - roc_obj$specificities,
  TPR = roc_obj$sensitivities
)

p <- ggplot(roc_df, aes(x = FPR, y = TPR)) +
  geom_line(color = "#8E44AD", size = 1.2) +   # 紫色区分
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  labs(title = "随机森林 ROC 曲线 (测试集)",
       subtitle = paste0("AUC = ", round(auc_val, 3), 
                         " (95% CI: ", round(ci_val[1], 3), " - ", round(ci_val[3], 3), ")"),
       x = "假阳性率 (1 - 特异度)",
       y = "真阳性率 (敏感度)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10)) +
  coord_equal() +
  xlim(0, 1) + ylim(0, 1)

ggsave("rf_roc_curve.pdf", plot = p, width = 6, height = 5, dpi = 300)
ggsave("rf_roc_curve.png", plot = p, width = 6, height = 5, dpi = 300)
print(p)

cat("\n✅ ROC 曲线已保存为：rf_roc_curve.pdf 和 rf_roc_curve.png\n")
cat("\n✅ 随机森林建模完成（手动10折CV + 完整指标 + ROC图）！\n")



# ==================== 加载包 ====================
library(pROC)
library(ggplot2)
library(dplyr)

# ==================== 1. 读取各模型预测结果 ====================
model_files <- list(
  "LR"   = "lr_test_predictions.csv",
  "SVM"  = "svm_test_predictions.csv",
  "DT"   = "dt_test_predictions.csv",
  "RF"   = "rf_test_predictions.csv"
)

# 模型中文全称（用于底部脚注）
model_full <- c("LR"="Logistic Regression", 
                "SVM"="Support Vector Machine", 
                "DT"="Decision Tree", 
                "RF"="Random Forest")

# 存储结果
roc_obj_list <- list()
auc_summary <- data.frame(Model = character(), 
                          AUC = numeric(), 
                          CI_low = numeric(), 
                          CI_high = numeric(),
                          stringsAsFactors = FALSE)

# ==================== 2. 计算每个模型的 AUC 和 95% CI ====================
for (name in names(model_files)) {
  dat <- read.csv(model_files[[name]])
  dat$true_label <- factor(dat$true_label, levels = c("No", "Yes"))
  
  # 计算ROC
  roc_obj <- roc(response = dat$true_label, predictor = dat$pred_prob, 
                 levels = c("No", "Yes"), quiet = TRUE)
  roc_obj_list[[name]] <- roc_obj
  
  # 计算AUC和95% CI（DeLong法，文献中最常用）
  ci_auc <- ci.auc(roc_obj, method = "delong")
  auc_summary <- rbind(auc_summary, data.frame(
    Model = name,
    AUC = round(ci_auc[2], 3),
    CI_low = round(ci_auc[1], 3),
    CI_high = round(ci_auc[3], 3)
  ))
}

# 构建图例标签：模型缩写 (AUC: 0.xxx, 95% CI: 0.xxx-0.xxx)
legend_labels <- paste0(
  auc_summary$Model, 
  " (AUC: ", sprintf("%.3f", auc_summary$AUC), 
  ", 95% CI: ", sprintf("%.3f", auc_summary$CI_low), 
  "-", sprintf("%.3f", auc_summary$CI_high), ")"
)

# ==================== 3. 将 ROC 数据转为绘图数据框 ====================
roc_df <- data.frame()
for (name in names(roc_obj_list)) {
  roc_obj <- roc_obj_list[[name]]
  # 提取坐标
  coords_data <- coords(roc_obj, "all", ret = c("fpr", "tpr"))
  roc_df <- rbind(roc_df, data.frame(
    fpr = coords_data$fpr,
    tpr = coords_data$tpr,
    model = name
  ))
}

# ==================== 4. 绘制文献风格 ROC 曲线 ====================
# 专业配色（色盲友好，打印清晰）
colors <- c(
  "LR"  = "#D55E00",   # 橙红（通常用于对照）
  "SVM" = "#0072B2",   # 深蓝
  "DT"  = "#009E73",   # 墨绿
  "RF"  = "#CC79A7"    # 粉紫
)

p <- ggplot(roc_df, aes(x = fpr, y = tpr, color = model)) +
  # ROC 曲线（粗实线）
  geom_line(size = 1.2) +
  
  # 对角线（随机猜测）
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", 
              color = "gray40", size = 0.8) +
  
  # 自定义颜色
  scale_color_manual(values = colors, labels = legend_labels) +
  
  # 坐标轴标签（中英文均可，这里保留英文更国际范）
  labs(
    x = "False Positive Rate (1 - Specificity)",
    y = "True Positive Rate (Sensitivity)",
    caption = paste0(
      "Note: ", 
      paste(names(model_full), model_full, sep = ": ", collapse = "; ")
    )
  ) +
  
  # 经典带灰色网格的主题（顶刊标配）
  theme_bw(base_size = 13) +
  
  # 精细调整
  theme(
    # === 灰色网格（您要求的） ===
    panel.grid.major = element_line(color = "gray85", size = 0.5),
    panel.grid.minor = element_line(color = "gray93", size = 0.3),
    
    # === 图例（置于底部，横向，含完整CI信息） ===
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 9.5, face = "bold"),
    legend.spacing.x = unit(0.2, "cm"),
    legend.box = "horizontal",
    legend.key.width = unit(1.2, "cm"),
    
    # === 坐标轴 ===
    axis.line = element_line(color = "black", size = 0.6),
    axis.text = element_text(color = "black", size = 11),
    axis.title = element_text(size = 12, face = "bold"),
    axis.ticks = element_line(color = "black", size = 0.6),
    
    # === 底部脚注（斜体，左对齐） ===
    plot.caption = element_text(hjust = 0, size = 9, face = "italic",
                                margin = margin(t = 10)),
    
    # === 边框（保留theme_bw的边框，更显正式） ===
    panel.border = element_rect(color = "black", fill = NA, size = 0.8),
    
    # 去除多余空白
    plot.margin = margin(10, 10, 15, 10)
  ) +
  
  # 保持1:1比例
  coord_equal() +
  xlim(0, 1) + ylim(0, 1)

# ==================== 5. 保存高分辨率图片 ====================
ggsave("ROC_文献风格_带网格_带CI脚注.pdf", plot = p, width = 8, height = 7, dpi = 300)
ggsave("ROC_文献风格_带网格_带CI脚注.png", plot = p, width = 8, height = 7, dpi = 300, bg = "white")

# 在 R 中显示
print(p)

# ==================== 6. 控制台输出汇总（方便您复制到论文表格） ====================
cat("\n========== AUC 汇总（含 95% CI） ==========\n")
print(auc_summary)

cat("\n✅ 图片已保存：")
cat("\n   - ROC_文献风格_带网格_带CI脚注.pdf")
cat("\n   - ROC_文献风格_带网格_带CI脚注.png")
cat("\n\n图例中已包含每个模型的 AUC 和 95% 置信区间。")
cat("\n底部脚注包含模型全称。\n")
